/* Author: Florian Seidler, Jonathan Zopes
 * Date: 25.11.2012
 * Issue: Misc Functions for Polymer Simulation
*/

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
#include <string.h>
#include "misc.h"

signed int randomsign(int n) {
    double uniformrv;
    uniformrv=((double)rand()/(double)RAND_MAX);
    if(n==0) {
        if(uniformrv>0.5) return -1;
        else return 1;
    }
    if(n<0) {
        return -1;
    }
    else {
        return 1;
    }
}

int generate_random_direction(int lastdirection, int dimension) {
    double uniformrv;
    int i;
    int direction=-lastdirection;
    while(direction==-lastdirection) {
        uniformrv=((double)rand()/(double)RAND_MAX);
        for(i=1; i<dimension+1;i++) {
            if(uniformrv<(double)(i)/(double)dimension) {
                direction=i;
                break;
            }
        }
        uniformrv=((double)rand()/(double)RAND_MAX);
        if(uniformrv<0.5) {
            direction*=-1;
        }
    }
    return direction;
}

double calc_walk_distance(int **walk,int walksize,int dimension) {
    int i;
    double distance=0;
    for(i=0;i<dimension;i++) {
        distance+=(walk[0][i]-walk[walksize-1][i])*(walk[0][i]-walk[walksize-1][i]);
    }
    return distance;
}

void copy_walk(int **fromwalk, int **towalk, int walksize, int dimension) {
    int i,j;
    for(i=0; i<walksize; i++) {
        for(j=0;j<dimension;j++) {
            towalk[i][j]=fromwalk[i][j];
        }
    }
}

void set_walk(int **setwalk, int walksize, int dimension, int c) {
    int i,j;
    for(i=0; i<walksize; i++) {
        for(j=0;j<dimension;j++) {
            setwalk[i][j]=c;
        }
    }
}

void set_array(int *array, int length, int c) {
    int i;
    for(i=0;i<length;i++) {
        array[i]=c;
    }
}

void print_walk_to_file(int **walk, int walksize, int dimension, char *path) {
    int i,j;
    FILE *fp;
    fp=fopen(path, "w");
     for(i=0;i<walksize;i++) {
        for(j=0;j<dimension;j++) {
            fprintf(fp,"%d ",walk[i][j]);
        }
        fprintf(fp,"\n");
    }
    fclose(fp);
}

int check_for_intersection(int **walk,int *buffer, int walksize, int dimension, int step) {
    int i,j,intersection;
    intersection=0;
    if(step>walksize) {
        printf("Wrong Stepsize. Bigger than walksize");
    }
    for(i=0;i<step;i++) {
        intersection=1;
        for(j=0;j<dimension;j++) {
            //Check if the new position to be taken was already taken before.
            if(walk[i][j]!=buffer[j]) {
                intersection=0;
            }
        }
        if (intersection==1) return 1;
    }
    return 0;
}

int simple_sampling_saw(int **walk,int *buffer, int walksize, int dimension) {
    int step,lastdirection,direction,intersection,rejects,i;
    step=1;
    rejects=0;
    lastdirection=0;
    intersection=0;
    while(step<walksize) {
        if(step==1) {
            /*Reset walks and afterwards
             *generate a  random direction with no lastdirection
             *since dimension+1 is no possible lastdirection.
            */
            set_walk(walk,walksize,dimension,0);
            direction=generate_random_direction(dimension+1,dimension);
            walk[1][abs(direction)-1]=walk[0][abs(direction)-1]+direction/abs(direction);
            lastdirection=direction;
            step++;
        }
        else {
            direction=generate_random_direction(lastdirection,dimension);
            for(i=0;i<dimension;i++) {
                if(i+1==abs(direction)) {
                    buffer[i]=walk[step-1][i]+direction/abs(direction);
                }
                else buffer[i]=walk[step-1][i];
            }
            intersection=check_for_intersection(walk,buffer,walksize,dimension,step);
            if(intersection==0) {
                for(i=0;i<dimension;i++) {
                   walk[step][i]=buffer[i];
                }
                step++;
                lastdirection=direction;
            }
            else {
                step=1;
                rejects++;
            }
        }
    }
    return rejects;
}

double biased_sampling_saw(int **walk, int *buffer, int *freedirection ,int walksize, int dimension) {
    int step,lastdirection,direction,intersection,rejects,freedirections,index,i,j;
    int rejected;
    double weight;
    step=1;
    rejects=0;
    lastdirection=0;
    intersection=0;
    while(step<walksize) {
        freedirections=0;
        set_array(freedirection,2*dimension,0);
        if(step==1) {
            weight=1;
            set_walk(walk,walksize,dimension,0);
            direction=generate_random_direction(dimension+1,dimension);
            walk[1][abs(direction)-1]=walk[0][abs(direction)-1]+direction/abs(direction);
            lastdirection=direction;
            step++;
        }else {
            //search for free directions
            freedirections=0;
            set_array(freedirection,2*dimension,0);
            for(i=0;i<dimension;i++) {
                for(j=0;j<dimension;j++) {
                    if(j==i) buffer[j]=walk[step-1][j]+1;
                    else buffer[j]=walk[step-1][j];
                }
                intersection=check_for_intersection(walk,buffer,walksize,dimension,step);
                if(intersection==0) {
                    freedirection[i]=1;
                    freedirections++;
                }
            }
            for(i=dimension;i<2*dimension;i++) {
                for(j=0;j<dimension;j++) {
                    if(j==i-dimension) buffer[j]=walk[step-1][j]-1;
                    else buffer[j]=walk[step-1][j];
                }
                intersection=check_for_intersection(walk,buffer,walksize,dimension,step);
                if(intersection==0) {
                    freedirection[i]=1;
                    freedirections++;
                }
            }
            if(freedirections!=0) {
                direction=0;
                while(direction==0) {
                    direction=generate_random_direction(dimension+10,dimension);
                    if(direction<0) {
                        index=abs(direction)-1+dimension;
                    }
                    else {
                        index=direction-1;
                    }
                    if(freedirection[index]==0) {
                        direction=0;
                    }
                }
                for(i=0;i<dimension;i++) {
                    if(i==abs(direction)-1) walk[step][i]=walk[step-1][i]+direction/abs(direction);
                    else walk[step][i]=walk[step-1][i];
                }
                step++;
                weight=weight*(double)freedirections/(double)(2*dimension-1);
            }else{
                step=1;
                rejects++;
            }
        }
    }
    return weight;
}

double simple_sampling_average_distance(int **walk, int *buffer, int walksize, int dimension, int numberofwalks) {
    int i;
    double average_distance=0;
    for(i=0;i<numberofwalks;i++) {
        simple_sampling_saw(walk,buffer,walksize,dimension);
        average_distance+=calc_walk_distance(walk,walksize,dimension);
    }
    return average_distance/(int)numberofwalks;
}





